set more off

**betai and zstar are unknown so we need to simulate them given z

local J `1'
local sigma `2'




*Create dataset with an observation for each (beneficiary, choice)

reshape long x_ z_ epsilon_ , i(ID c alpha beta) j(altID)

rename x_ x
rename z_ z
rename epsilon_ epsilon

*Solve for reservation values as a function of costs
*Kim et. al. (2010)
* c / (beta*sigma) = normpdf(r/sigma)-(r/sigma)*(1-normcdf(r/sigma))
* sigma is the standard deviation of z

do programs/solvefor_r.do

capture program drop my_fn
program my_fn
    syntax, fn(varname) fp(varname) x(varname) a(varname)
    replace `fn' = normalden(`x') - `x' + `x' * normal(`x') - `a'
    replace `fp' = normal(`x') - 1
end

*a is c/(beta*sigma)

gen a = c/(beta*`sigma')
gen x0 = -a
my_nr rdivsigma, a(a) x0(x0) max(1000) disp(50) f(my_fn)
egen meanz = mean(z)
gen r = beta*rdivsigma*`sigma'+x*alpha+epsilon+beta*meanz
drop a x0 rdivsigma meanz

*Create choice and open numbers
*Order choices by free utility
*Keep searching if:
*max_{k \in set of searched goods}  alpha*x_{k} + beta*z_{k} + epsilon_{k}   < alpha*x_{j+1} + beta*r + epsilon_{j+1}

*Loop through goods
*One variable is current best
*If current best exceeds benchmark for next good, stop and choose current best
*If current best is less than benchmark, keep going

*At the end of loop, need to tag open boxes and chosen good

gen FU = alpha*x+epsilon
gen fullU = FU+beta*z
*Note, r is now defined in utility space
gen reserve = r
*gen reserve = FU+beta*r

gsort ID -FU
by ID: gen newID = _n

forvalues x = 1/`J' {
by ID: gen FU_`x' = FU[`x']
by ID: gen U_`x' = fullU[`x']
by ID: gen reserve_`x' = reserve[`x']
}

gen open = 1
gen Ucurr = U_1
local Jminus = `J'-1
forvalues x = 1/`Jminus' {
local z = `x'+1
replace open = `z' if Ucurr < reserve_`z'
*If so, update utility if its bigger
replace Ucurr = U_`z' if U_`z' > Ucurr & open == `z'
}

*Now open is correct
gen choice = 0
forvalues x = 1/`J' {
replace choice = `x' if Ucurr == U_`x' & `x' <= open
}


rename open numopen
gen open = newID <= numopen
gen chosen = (choice == newID)

bys ID: egen numchosen = sum(chosen)
egen maxchosen = max(numchosen)
local maxchosen = maxchosen
while(`maxchosen' != 1) {
drop numchosen maxchosen
replace tempU = tempU+.00001*uniform()
drop chosen
gen chosen = (tempU == maxU)
bys ID: egen numchosen = sum(chosen)
egen maxchosen = max(numchosen)
local maxchosen = maxchosen
}
drop numchosen maxchosen

drop U_* reserve_* FU_*
